Se requieren cargar las siguientes paqueterĆas para seguir los ejemplos.
library(easypackages)
libraries("tidyverse","fpp3","plotly", "patchwork")Los modelos que veremos aquà tienen como idea principal encontrar relaciones lineales entre la serie que queremos pronosticar, \(y\), con una o mÔs series distintas, x. En otras palabras, pronosticaremos los valores futuros de una serie, a partir de los cambios en otra serie que la afecte.
Es muy comĆŗn querer predecir de esta forma. Por ejemplo, una tienda de helados podrĆa encontrar una relación entre sus ventas (\(y\)) y la temperatura (\(x_1\)). O las ventas de Nike, a partir de cuĆ”nto gastan en publicidad y mercadotecnia.
En la literatura podemos encontrar muchos nombres para las variables \(y\) ^ \(x\). P. ej.
| \(y\) (var. de pronóstico) | \(x\) (vars. predictoras) |
|---|---|
| Var. dependiente | Vars. independientes |
| Explicada | Explicativas |
| Regresada | Regresoras |
| Respuesta | EstĆmulo |
| Resultado | Covariante |
| Controlada | De control |
El caso mĆ”s sencillo serĆa un modelo de regresión lineal simple, de la forma:
\[ y_t = \beta_0 + \beta_1 x_t + \varepsilon_t \]
donde
\(\beta_0\) es conocido como el intercepto y representa el valor predicho cuando \(x = 0\).
\(\beta_1\) es la pendiente de la recta. Nos indica el cambio promedio en \(y\), ante un cambio en una unidad de \(x\).
El tĆ©rmino de error, \(\varepsilon_t\) se asume aleatorio y decimos que captura los cambios debido a todas las otras variables que pudieran llegar a afectar a \(y_t\), que no estĆ”n explĆcitamente especificadas en el modelo.
La recta resultante estƔ dada entonces por \(\beta_0 + \beta_1 x_t\), y la diferencia que existe en los puntos reales y Ʃsta es \(\varepsilon_t\).
Como primer ejemplo, veamos las tasas de crecimiento del gasto de consumo, \(y\), y su relación con el ingreso personal disponible, \(x\).
La grƔfica de tiempo de ambas series:
us_change %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Consumption, colour = "Consumo")) +
geom_line(aes(y = Income, colour = "Ingreso")) +
ylab("cambio %") + xlab("AƱo") +
guides(colour=guide_legend(title="Series")) +
theme(legend.position = "top")Un diagrama de dispersión entre ambas series, para ver una posible correlación.
us_change %>%
ggplot(aes(x=Income, y=Consumption)) +
ylab("Consumo (cambio % trimestral)") +
xlab("Ingreso (cambio % trimestral)") +
geom_point() +
geom_smooth(method="lm", se=FALSE)La lĆnea azul en el grĆ”fico sigue la ecuación que describe la regresión lineal que tiene por variable dependiente al consumo y como independiente al ingreso.
El tĆ©rmino de regresión fue acuƱado por primera vez por Francis Galton en 1886. Ćl estaba estudiando la relación que existe entre la estatura de los hijos y la estatura de los padres.
Lo que encontró fue lo siguiente, en resumen:
Los padres mĆ”s altos, tendĆan a tener hijos mĆ”s altos, mientras que los padres bajos tendĆan a tener hijos bajos.
En promedio, los hijos de padres altos no logran ser mƔs altos que ellos. Similarmente, los hijos de padres bajos, en promedio son mƔs altos que sus papƔs.
AsĆ, Galton decĆa que habĆa una tendencia a regresar a la estatura promedio.
De no cumplirse la regresión de Galton, serĆa comĆŗn tener gente de la estatura de un Hobbit, y tambiĆ©n de la estatura de un gigante.
Entonces, el anÔlisis de regresión en tiempos modernos trata sobre la relación de la dependencia entre una variable \(y\), respecto de una o mÔs variables exógenas (regresoras \(x\)) para predecir el valor promedio de la variable dependiente.
āUna relación estadĆstica, por mĆ”s fuerte y sugerente que sea, nunca podrĆ” establecer una conexión causal: nuestras ideas de causalidad deben provenir de estadĆsticas externas y, en Ćŗltimo tĆ©rmino, de una u otra teorĆaā (Kendall & Stuart, 1961)
Regresión \(\neq\) Causalidad
Correlación \(\neq\) Causalidad
Se puede hablar de linealidad en dos sentidos:
Para un modelo de regresión lineal, solo nos interesa que sea lineal en los parÔmetros.
Todas estas funciones son lineales en los parÔmetros y pueden ser estimadas mediante un modelo de regresión lineal
AsĆ, un modelo de regresión lineal puede generar una recta, o una variedad de curvas, dependiendo la forma funcional que se elija.
\[ H_0: \beta_0 = 0 \\ H_1: \beta_0 \neq 0 \] \[ H_0: \beta_1 = 0 \\ H_1: \beta_1 \neq 0 \]
\[ \hat{y}_t = \hat{\beta}_0 + \hat{\beta}_1x_t + \hat{\varepsilon}_t \\ \hat{y}_t = 0.54454 + 0.27183 x_t + \hat{\varepsilon}_t \]
\[ y_{consumo} = \beta_0 + \beta_1 x_{income} + \beta_2 x_{production} + \beta_3 x_{savings} + \beta_4 x_{unemployment} \]
us_changeus_change %>%
as_tibble() %>%
select(-Quarter) %>%
GGally::ggpairs()Una correlación, por mÔs fuerte que sea entre dos variables, no puede implicar por sà misma causalidad.
us_change %>%
pivot_longer(cols = -Quarter) %>%
ggplot(aes(x = Quarter, y = value, color = name)) +
geom_line() +
facet_wrap(~ name, scales = "free_y") +
theme(legend.position = "none")Graficando nuestra variable de pronóstico (Consumo) vs. cada una de las variables predictoras:
us_change %>%
pivot_longer(cols = -c(Quarter, Consumption)) %>%
ggplot(aes(x = Quarter, y = value, color = name)) +
geom_line() +
geom_line(aes(y = Consumption), color = "black") +
facet_wrap(~ name, scales = "free_y") +
theme(legend.position = "none")Realizamos un primer modelo, donde utilizaremos de variable predictora al ingreso disponible, para pronosticar el consumo.
\[ y_{t,Consumo} = \beta_0 + \beta_{1}x_{t,Ingreso} + \varepsilon_t \]
fit1 <- us_change %>%
model(reg_lin_simple = TSLM(Consumption ~ Income)
)
fit1 %>% report()Series: Consumption
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-2.58236 -0.27777 0.01862 0.32330 1.42229
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.54454 0.05403 10.079 < 2e-16 ***
Income 0.27183 0.04673 5.817 2.4e-08 ***
---
Signif. codes:
0 ā***ā 0.001 ā**ā 0.01 ā*ā 0.05 ā.ā 0.1 ā ā 1
Residual standard error: 0.5905 on 196 degrees of freedom
Multiple R-squared: 0.1472, Adjusted R-squared: 0.1429
F-statistic: 33.84 on 1 and 196 DF, p-value: 2.4022e-08
Prueba de significancia individual para cada uno de los parƔmetros (\(\beta_i\))
\[ H_0: \beta_i = 0 \]
Prueba de significancia conjunta (para revisar si nuestro modelo sirve o no)
\[ H_0: \beta_1 = \beta_2 = \beta_3 = \ldots = 0 \]
augment(fit1) %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Consumption, color = "Datos")) +
geom_line(aes(y = .fitted, color = "Fitted"))+
xlab("AƱo") + ylab(NULL) +
ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
guides(color = guide_legend(title = NULL))El modelo no parece capturar adecuadamente la variación de los datos reales.
augment(fit1) %>%
ggplot(aes(x = Consumption, y = .fitted)) +
geom_point() +
ylab("Fitted (valores ajustados)") +
xlab("Datos (reales históricos)") +
ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
geom_abline(intercept = 0, slope = 1)fit1 %>%
gg_tsresiduals()augment(fit1) %>%
features(.resid, ljung_box, lag= 10, dof = 2)Es evidente que este modelo se puede mejorar. Probemos incluyendo las otras predictoras.
PodrĆamos proponer ahora un modelo en donde el consumo sea una función del ingreso, la producción, los ahorros y el desempleo:
\[ y_{t,Consumo} = \beta_0 + \beta_{1}x_{t,Ingreso} + \beta_2x_{t,Producción} + \beta_3x_{t,Ahorros} + \beta_4x_{t,Desempleo} + \varepsilon_t \]
fit2 <- us_change %>%
model(
reg_lin_multiple = TSLM(Consumption ~ Income + Production + Savings + Unemployment)
)
report(fit2)Series: Consumption
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-0.90555 -0.15821 -0.03608 0.13618 1.15471
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.253105 0.034470 7.343 5.71e-12 ***
Income 0.740583 0.040115 18.461 < 2e-16 ***
Production 0.047173 0.023142 2.038 0.0429 *
Savings -0.052890 0.002924 -18.088 < 2e-16 ***
Unemployment -0.174685 0.095511 -1.829 0.0689 .
---
Signif. codes:
0 ā***ā 0.001 ā**ā 0.01 ā*ā 0.05 ā.ā 0.1 ā ā 1
Residual standard error: 0.3102 on 193 degrees of freedom
Multiple R-squared: 0.7683, Adjusted R-squared: 0.7635
F-statistic: 160 on 4 and 193 DF, p-value: < 2.22e-16
augment(fit2) %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Consumption, color = "Datos")) +
geom_line(aes(y = .fitted, color = "Fitted"))+
xlab("AƱo") + ylab(NULL) +
ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
guides(color = guide_legend(title = NULL))Este modelo parece capturar mÔs variación de los datos históricos.
augment(fit2) %>%
ggplot(aes(x = Consumption, y = .fitted)) +
geom_point() +
ylab("Fitted (valores ajustados)") +
xlab("Datos (reales históricos)") +
ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
geom_abline(intercept = 0, slope = 1)fit2 %>%
gg_tsresiduals()augment(fit2) %>%
features(.resid, ljung_box, lag= 10, dof = 2)df <- left_join(us_change, residuals(fit2), by = "Quarter")
df %>%
select(-c(Consumption, .model)) %>%
pivot_longer(cols = c(Income:Unemployment)) %>%
ggplot(aes( x = value, y = .resid, color = name)) +
geom_point() + ylab("Residuales") + xlab("Predictoras") +
facet_wrap(~ name, scales = "free_x") +
theme(legend.position = "none")augment(fit2) %>%
ggplot(aes(x = .fitted, y = .resid)) +
geom_point() +
labs(x = "Ajustados", y = "Residuales")glance(fit2) %>%
select(adj_r_squared, AIC, AICc, BIC)fit3 <- us_change %>%
model(r1 = TSLM(Consumption ~ Income),
r2 = TSLM(Consumption ~ Income + Production),
r3 = TSLM(Consumption ~ Income + Production + Savings + Unemployment),
r4 = TSLM(Consumption ~ Income + Production + Savings),
r5 = TSLM(Consumption ~ Income + Savings + Unemployment),
r6 = TSLM(Consumption ~ Income + Production + Unemployment),
r7 = TSLM(Consumption ~ Income + Savings)
)
fit3 %>%
glance() %>%
select(.model, adj_r_squared, AIC, AICc, BIC)fit3 %>%
select(r3) %>%
report()Series: Consumption
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-0.90555 -0.15821 -0.03608 0.13618 1.15471
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.253105 0.034470 7.343 5.71e-12 ***
Income 0.740583 0.040115 18.461 < 2e-16 ***
Production 0.047173 0.023142 2.038 0.0429 *
Savings -0.052890 0.002924 -18.088 < 2e-16 ***
Unemployment -0.174685 0.095511 -1.829 0.0689 .
---
Signif. codes:
0 ā***ā 0.001 ā**ā 0.01 ā*ā 0.05 ā.ā 0.1 ā ā 1
Residual standard error: 0.3102 on 193 degrees of freedom
Multiple R-squared: 0.7683, Adjusted R-squared: 0.7635
F-statistic: 160 on 4 and 193 DF, p-value: < 2.22e-16
Selección del modelo a través de backwards stepwise regression con base en el \(AICc\) sugiere conservar el modelo que incluye todas las predictoras.
us_change %>%
model(i = TSLM(Consumption ~ Income + Production + Savings + Unemployment),
ii = TSLM(Consumption ~ Income + Production + Savings),
iii = TSLM(Consumption ~ Income + Production + Unemployment),
iv = TSLM(Consumption ~ Income + Savings + Unemployment),
v = TSLM(Consumption ~ Production + Savings + Unemployment)
) %>%
glance() %>%
select(.model, adj_r_squared, AIC, AICc, BIC)Con base en \(BIC\) parece que el mejor modelo es el que incluye todas las predictoras menos el desempleo.
us_change %>%
model(i = TSLM(Consumption ~ Income + Production + Savings + Unemployment),
ii = TSLM(Consumption ~ Income + Production + Savings),
iii = TSLM(Consumption ~ Income + Production),
iv = TSLM(Consumption ~ Income + Savings),
v = TSLM(Consumption ~ Production + Savings),
vi = TSLM(Consumption ~ Income + Savings + Unemployment),
vii = TSLM(Consumption ~ Production + Savings + Unemployment),
viii = TSLM(Consumption ~ Income + Production + Unemployment)
) %>%
glance() %>%
select(.model, adj_r_squared, AIC, AICc, BIC)us_change %>%
model(i = TSLM(Consumption ~ Income),
ii = TSLM(Consumption ~ Production),
iii = TSLM(Consumption ~ Savings),
iv = TSLM(Consumption ~ Unemployment)
) %>%
glance() %>%
select(.model, adj_r_squared, AIC, AICc, BIC)El mejor modelo hasta ahora es el que toma de predictora a la producción.
us_change %>%
model(
o = TSLM(Consumption ~ Production),
i = TSLM(Consumption ~ Production + Income),
ii = TSLM(Consumption ~ Production + Savings),
iii = TSLM(Consumption ~ Production + Unemployment)
) %>%
glance() %>%
select(.model, adj_r_squared, AIC, AICc, BIC)El mejor hasta ahora es el modelo con la producción y el ingreso.
us_change %>%
model(
i = TSLM(Consumption ~ Production + Income),
ii = TSLM(Consumption ~ Production + Income + Savings),
iii = TSLM(Consumption ~ Production + Income + Unemployment)
) %>%
glance() %>%
select(.model, adj_r_squared, AIC, AICc, BIC)El modelo mejora agregando los ahorros tambiƩn.
us_change %>%
model(
i = TSLM(Consumption ~ Production + Income + Savings),
ii = TSLM(Consumption ~ Production + Income + Savings + Unemployment)
) %>%
glance() %>%
select(.model, adj_r_squared, AIC, AICc, BIC)Utilizando la adj_r_squared o el AICc, concluimos lo mismo que antes: el modelo que incluye a todas las variables es el mejor. ### {-}
En estos pronósticos solo se utiliza información disponible hasta el último dato del histórico. A estos pronósticos se les considera como pronósticos reales. Aquà las predictoras se deben pronosticar antes de poder producir el pronóstico de la variable de interés.
Con estos pronósticos se utiliza información real disponible de las predictoras. Estos pronósticos ya no son reales (en el sentido estricto). La variable a pronosticar (\(y\)) sigue siendo desconocida.
fit_escenarios <- us_change %>%
model(lineal = TSLM(Consumption ~ Income + Savings + Unemployment))
# Necesitamos agregar nuevos datos de las predictoras
escenarios <- scenarios(
optimista = new_data(us_change, 4) %>%
mutate(Income = c(0,0.2,0,1.2), Savings = c(0,-0.1,0.5,-1), Unemployment = -0.1),
pesimista = new_data(us_change, 4) %>%
mutate(Income = -1, Savings = -0.5, Unemployment = 0.1),
names_to = "Escenario"
)
fc_escenarios <- fit_escenarios %>%
forecast(new_data = escenarios)
us_change %>%
autoplot(Consumption) +
autolayer(fc_escenarios)Vamos a generar un pronóstico ex-ante del consumo. Para esto, necesitamos primero producir pronósticos de las predictoras:
# mod_predictoras <- function(predictora, horizonte = 4) {
# us_change %>%
# model(predictora = ARIMA(as.formula(predictora)) %>%
# forecast(h = horizonte)
# }
#
# mod_predictoras(predictora = Income)
ingreso <- us_change %>%
model(ETS = ETS(Income),
ARIMA = ARIMA(Income)
) %>%
forecast(h = 4)
ingreso %>%
autoplot(us_change, level = NULL)
produccion <- us_change %>%
model(
ETS = ETS(Production),
ARIMA = ARIMA(Production)
) %>%
forecast(h = 4)
produccion %>%
autoplot(us_change, level = NULL)
ahorro <- us_change %>%
model(
ETS = ETS(Savings),
ARIMA = ARIMA(Savings)
) %>%
forecast(h = 4)
ahorro %>%
autoplot(us_change, level = NULL)
desempleo <- us_change %>%
model(
ETS = ETS(Unemployment),
ARIMA = ARIMA(Unemployment)
) %>%
forecast(h = 4)
desempleo %>%
autoplot(us_change, level = NULL)Teniendo ya los pronósticos de cada predictora, podemos proceder a generar el pronóstico del consumo.
fit <- us_change %>%
model(
`Regresión lineal múltiple` = TSLM(Consumption ~ Income + Production + Savings + Unemployment)
)
datos_futuros <- new_data(us_change,4) %>%
mutate(Income = ingreso %>% filter(.model == "ARIMA") %>% pull(.mean),
Savings = ahorro %>% filter(.model == "ARIMA") %>% pull(.mean),
Unemployment = desempleo %>% filter(.model == "ARIMA") %>% pull(.mean),
Production = produccion %>% filter(.model == "ARIMA") %>% pull(.mean))
datos_futuros
fc <- forecast(fit, datos_futuros)
fc %>%
autoplot(us_change)fc %>%
autoplot(us_change %>% filter_index("2016 Q1" ~ .))recent_production <- aus_production %>%
filter(year(Quarter) >= 1992)
recent_production
recent_production %>%
autoplot(Beer) +
labs(x = "AƱo", y = "Megalitros",
title = "Producción de cerveza trimestral en Australia")Existen varias predictoras que pueden ser útiles en el anÔlisis de regresión.
trend()\[ y_t = \beta_0 + \beta_1t + \varepsilon_t \]
Las variables dummy (tambiĆ©n llamadas dicotómicas, binarias, ā¦) toman solo dos valores: 1 o 0 (verdadero/falso).
Para agregar variables dummy estacionales, basta con escribir season().
recent_production %>%
select(Quarter,Beer) %>%
mutate(tendencia = seq_along(recent_production$Quarter),
q2 = if_else(quarter(Quarter)==2,1,0),
q3 = if_else(quarter(Quarter)==3,1,0),
q4 = if_else(quarter(Quarter)==4,1,0)
) %>%
model(TSLM(Beer ~ tendencia + q2 + q3 + q4)) %>%
report()Series: Beer
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-42.9029 -7.5995 -0.4594 7.9908 21.7895
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 441.80044 3.73353 118.333 < 2e-16 ***
tendencia -0.34027 0.06657 -5.111 2.73e-06 ***
q2 -34.65973 3.96832 -8.734 9.10e-13 ***
q3 -17.82164 4.02249 -4.430 3.45e-05 ***
q4 72.79641 4.02305 18.095 < 2e-16 ***
---
Signif. codes:
0 ā***ā 0.001 ā**ā 0.01 ā*ā 0.05 ā.ā 0.1 ā ā 1
Residual standard error: 12.23 on 69 degrees of freedom
Multiple R-squared: 0.9243, Adjusted R-squared: 0.9199
F-statistic: 210.7 on 4 and 69 DF, p-value: < 2.22e-16
fit_beer <- recent_production %>%
model(TSLM(Beer ~ trend() + season()))
report(fit_beer)Series: Beer
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-42.9029 -7.5995 -0.4594 7.9908 21.7895
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 441.80044 3.73353 118.333 < 2e-16
trend() -0.34027 0.06657 -5.111 2.73e-06
season()year2 -34.65973 3.96832 -8.734 9.10e-13
season()year3 -17.82164 4.02249 -4.430 3.45e-05
season()year4 72.79641 4.02305 18.095 < 2e-16
(Intercept) ***
trend() ***
season()year2 ***
season()year3 ***
season()year4 ***
---
Signif. codes:
0 ā***ā 0.001 ā**ā 0.01 ā*ā 0.05 ā.ā 0.1 ā ā 1
Residual standard error: 12.23 on 69 degrees of freedom
Multiple R-squared: 0.9243, Adjusted R-squared: 0.9199
F-statistic: 210.7 on 4 and 69 DF, p-value: < 2.22e-16
p <- augment(fit_beer) %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Beer, color = "Datos")) +
geom_line(aes(y = .fitted, color = "Fitted")) +
labs(x = "AƱo", y = "Megalitros")
ggplotly(p)Al ver la grÔfica, vemos que existe un trimestre muy por debajo del resto. Podemos agregar una variable de intervención para ese periodo.
Capturan el efecto de un solo periodo.
cerveza <- recent_production %>%
select(Quarter, Beer) %>%
mutate(
q2_94 = if_else(Quarter == yearquarter("1994 Q2"),1,0),
q4_04 = if_else(Quarter == yearquarter("2004 Q4"),1,0)
)
cervezaAjustamos un modelo, corrigiendo por ese periodo outlier (1994 Q2).
fit_beer <- cerveza %>%
model(TSLM(Beer ~ trend() + season() + q2_94 + q4_04))
report(fit_beer)Series: Beer
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-27.0379 -5.8811 -0.6895 7.7209 21.7895
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 441.84123 3.32189 133.009 < 2e-16
trend() -0.34137 0.05975 -5.713 2.78e-07
season()year2 -33.39369 3.55788 -9.386 7.83e-14
season()year3 -17.82164 3.55460 -5.014 4.16e-06
season()year4 75.32030 3.60790 20.876 < 2e-16
q2_94 -24.03384 11.24265 -2.138 0.036190
q4_04 -45.41027 11.15547 -4.071 0.000126
(Intercept) ***
trend() ***
season()year2 ***
season()year3 ***
season()year4 ***
q2_94 *
q4_04 ***
---
Signif. codes:
0 ā***ā 0.001 ā**ā 0.01 ā*ā 0.05 ā.ā 0.1 ā ā 1
Residual standard error: 10.81 on 67 degrees of freedom
Multiple R-squared: 0.9426, Adjusted R-squared: 0.9375
F-statistic: 183.4 on 6 and 67 DF, p-value: < 2.22e-16
p <- augment(fit_beer) %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Beer, color = "Datos")) +
geom_line(aes(y = .fitted, color = "Fitted")) +
labs(x = "AƱo", y = "Megalitros")
ggplotly(p)Capturan el efecto a partir de cierto periodo.
cerveza <- cerveza %>%
mutate(d2000 = if_else(year(Quarter)>=2000,1,0))
cervezafit_beer <- cerveza %>%
model(TSLM(Beer ~ trend() + season() + d2000))
report(fit_beer)Series: Beer
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-43.0235 -7.5945 -0.4539 7.7754 21.4829
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 441.9154 3.8644 114.354 < 2e-16 ***
trend() -0.3548 0.1308 -2.712 0.00846 **
season()year2 -34.6452 3.9985 -8.665 1.36e-12 ***
season()year3 -17.8046 4.0536 -4.392 4.02e-05 ***
season()year4 72.8279 4.0594 17.941 < 2e-16 ***
d2000 0.7282 5.6402 0.129 0.89766
---
Signif. codes:
0 ā***ā 0.001 ā**ā 0.01 ā*ā 0.05 ā.ā 0.1 ā ā 1
Residual standard error: 12.32 on 68 degrees of freedom
Multiple R-squared: 0.9243, Adjusted R-squared: 0.9188
F-statistic: 166.1 on 5 and 68 DF, p-value: < 2.22e-16
p <- augment(fit_beer) %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Beer, color = "Datos")) +
geom_line(aes(y = .fitted, color = "Fitted")) +
labs(x = "AƱo", y = "Megalitros")
ggplotly(p)cerveza %>%
gg_tsdisplay(Beer %>% difference(4), plot_type = "partial")tibble(mujer = c(1,0,1, 1),
hombre = c(0,1,0, 0),
nombre = c("Andrea","Juan","SofĆa","Fer"))Cuando se crean varaibles dummy, la cantidad de dummies a generar es una menos que el total de categorĆas. Si son dos categorĆas, solo necesitamos una dummy. Si son tres, necesitamos 2, etc.
\[ \log{y_t} = \beta_0 + \beta_1 \log{x_t} \]
La interpretación de los coeficientes (\(\beta_s\)) es como elasticidades (cambios porcentuales).
\[ y_t = \beta_0 + \beta_1 \log{x_t} \]
\[ \log{y_t} = \beta_0 + \beta_1 x_t \]
Aquà la regresión se calcula por partes en el tiempo, para capturar distintas tendencias.
boston_men <- boston_marathon %>%
filter(Event == "Men's open division") %>%
mutate(Minutes = as.numeric(Time)/60)
p <- boston_men %>%
autoplot(Minutes) +
ggtitle("Tiempos ganadores del maratón de Boston, categorĆa abierta de hombres")
ggplotly(p)fit_boston <- boston_men %>%
model(
lineal = TSLM(Minutes ~ trend()),
)
fc_boston <- fit_boston %>% forecast(h = 10)
boston_men %>%
autoplot(Minutes) +
geom_line(aes(y = .fitted, color = .model), data = fitted(fit_boston)) +
autolayer(fc_boston, alpha = 0.5, level = 95) +
ggtitle("Maratón de Boston, cat. abierta de hombres")Vamos a modelar los tiempos con una regresión por partes.
fit_boston <- boston_men %>%
model(
lineal = TSLM(Minutes ~ trend()),
exponencial = TSLM(log(Minutes) ~ trend()),
`Reg. por partes` = TSLM(Minutes ~ trend(knots = c(1940,1980)))
)
fc_boston <- fit_boston %>% forecast(h = 10)
boston_men %>%
autoplot(Minutes) +
geom_line(aes(y = .fitted, color = .model), data = fitted(fit_boston)) +
autolayer(fc_boston, alpha = 0.5, level = 95) +
ggtitle("Maratón de Boston, cat. abierta de hombres")Intentamos ahora poniendo los cambios en la tendencia en otros periodos:
fit_boston <- boston_men %>%
model(
lineal = TSLM(Minutes ~ trend()),
exponencial = TSLM(log(Minutes) ~ trend()),
`Reg. por partes` = TSLM(Minutes ~ trend(knots = c(1940,1980))),
`Reg. por partes2` = TSLM(Minutes ~ trend(knots = c(1975))),
`Reg. por partes3` = TSLM(Minutes ~ trend(knots = c(1915, 1950, 1988))),
`Reg. por partes4` = TSLM(Minutes ~ trend(knots = c(1915, 1927, 1950, 1985)))
)
fc_boston <- fit_boston %>% forecast(h = 10)
boston_men %>%
autoplot(Minutes) +
geom_line(aes(y = .fitted, color = .model), data = fitted(fit_boston)) +
autolayer(fc_boston, alpha = 0.5, level = 95) +
ggtitle("Maratón de Boston, cat. abierta de hombres")accuracy(fit_boston) %>%
arrange(MAPE)Probamos transformando la serie con Box-Cox:
lambda <- boston_men %>%
features(Minutes, guerrero) %>%
pull(lambda_guerrero)
p1 <- boston_men %>%
autoplot(Minutes)
p2 <- boston_men %>%
autoplot(box_cox(Minutes, lambda = lambda))
p3 <- boston_men %>%
autoplot(log(Minutes))
p1/p2/p3